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AN  ANALYTICAL  AND  NUMERICAL  INVESTIGATION 
ON  FAILURE  WAVES 


Abstract 

When  a  glass  specimen  is  shocked  near  but  below  the  apparent  Hugoniot  elastic 
limit,  it  may  undergo  elastic  deformations  at  the  shock  wave  front,  and  fail 
catastrophically  at  a  later  time.  Because  such  a  phenomenon  looks  different  from  usual 
inelastic  waves,  it  has  been  interpreted  as  a  failure  wave.  Based  on  the  observation  that 
the  failure  wave  propagates  with  a  degraded  speed  at  some  distance  behind  the  elastic 
shock  wave  front,  it  has  been  proposed  that  the  progressive  percolation  of  microfissures 
into  the  material  bulk  governs  the  failure  wave  mechanisms.  Stress  concentration  due  to 
the  defects  and  transient  loading  conditions  on  the  impact  surface  has  been  assumed  to  be 
the  origin  for  initiating  the  evolution  of  heterogeneous  microdamage.  To  provide  a 
general  framework  for  modeling  the  impact  response  of  certain  materials,  a  three- 
dimensional  continuum  damage  model  is  developed  in  this  report,  based  on  the 
assumption  of  shear-induced  dilatancy.  The  essential  component  of  the  proposed  model  is 
that  the  deviatoric  potential  energy  in  the  intact  material  is  converted  into  the  volumetric 
potential  energy  in  the  comminuted  and  dilated  material  during  the  time-dependent 
failure  evolution  process.  The  progressive  percolation  of  microfissures  is  described  by  a 
nonlinear  diffusion  equation  throughout  the  continuum  body.  The  diffusion  equation  and 
wave  equation  are  then  solved  via  a  staggered  manner  in  a  single  computational  domain. 
Numerical  solutions  are  presented  and  verified  with  experimental  data  available.  It 
appears  that  the  essential  feature  of  the  failure  wave  phenomenon,  as  observed  in  shock 
experiments  on  glasses,  can  be  predicted  by  the  proposed  approach.  A  very  recent 
experimental  study  indicates  that  the  dynamic  failure  evolution  in  mortar  is  a  rather 
gradual  process,  in  contrast  to  the  well-defined  failure  wave  front  as  observed  in  shocked 
glasses.  However,  the  lack  of  a  consistent  set  of  experimental  data  does  not  warrant  the 
revision  of  the  proposed  constitutive  model  to  predict  the  failure  wave  in  mortar. 
Although  a  considerable  progress  has  been  made  in  understanding  the  failure  wave 
phenomenon  since  it  was  reported  about  10  years  ago,  a  combined  experimental, 
analytical  and  computational  effort  is  still  needed  to  fully  understand  the  physics 
involved  in  the  dynamic  failure  behaviors  of  brittle  solids  under  various  impact  loads. 
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1.  Introduction 


When  a  glass  specimen  is  shocked  near  but  below  the  apparent  Hugoniot  elastic 
limit  (HEL),  it  may  undergo  elastic  deformations  at  the  shock  wave  front,  and  fail 
catastrophically  at  a  later  time.  Because  such  a  phenomenon  looks  different  from  usual 
inelastic  waves,  it  has  been  interpreted  as  a  failure  wave.  Since  Brar  et  al.  (1991)  and 
Kanel  et  al.  (1991)  reported  the  formation  and  propagation  of  failure  waves  in  a  series  of 
impact  experiments  with  glass  plates  and  bars,  continued  efforts  have  been  made  to 
explore  this  interesting  physical  phenomenon  (Bless  and  Brar,  1994;  Bourne  et  al.,  1995; 
Chen  and  Xin,  1999;  Clifton,  1993;  Espinosa  et  alM  1997a  and  b;  Feng,  2000;  Grady, 
1995a  and  b;  Grote  et  al.,  2000;  Raiser  and  Clifton,  1994;  Raiser  et  al.,  1994;  Rosenberg 
et  al.,  1996;  among  others).  However,  no  consensus  can  be  made  at  the  moment  on  the 
physics  behind  the  failure  wave  phenomenon  in  brittle  solids  under  impact  loading.  To 
set  a  stage  for  presenting  the  proposed  damage  model,  the  difference  between  a  usual 
inelastic  wave  and  the  failure  wave  is  first  analyzed,  based  on  the  literature  available,  as 
below. 

The  shock  response  of  glasses  beyond  the  HEL  often  displays  a  well-defined  two- 
wave  structure  that  has  been  interpreted  as  an  inelastic  response.  Although  the  prompt 
microcracking  across  the  compression  wave  front  occurs  if  the  glass  specimen  is  shocked 
at  or  above  the  HEL,  the  macroscopic  two-wave  structure  has  suggested  that  the 
compressive  failure  might  be  suppressed  by  confining  pressure  at  a  shocked  state  beyond 
the  HEL.  The  failure  wave  phenomenon  of  glasses,  which  occurs  when  the  compressive 
shock  stress  is  near  but  below  the  HEL,  suggests  that  the  HEL  may  not  be  an  elastic  limit, 
but  rather,  may  be  a  transition  in  failure  mechanisms.  A  possible  transition  is  the  one 
from  a  delayed  kinetic-controlled  failure  process  below  the  HEL  to  a  prompt  stress- 
controlled  failure  process  above  the  HEL  (Grady,  1995b).  Another  possibility  is  that  the 
HEL  may  represent  the  stress  level  above  which  bulk  glass  undergoes  permanent 
densification  (Espinosa  et  al.,  1997a).  In  other  words,  the  density  of  glass  rubble  due  to 
the  failure  below  the  HEL  is  different  from  that  above  the  HEL.  The  signature  feature 
that  separates  the  failure  wave  from  the  usual  inelastic  shock  wave  in  brittle  solids  is  that 
only  the  lateral  stress  is  changed  significantly  across  the  failure  front,  while  the 
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longitudinal  stress  is  almost  constant  during  the  failure  propagation  (Brar  et  al.,  1991; 
Kanel  et  al.,  1991;  Bourne  et  al.,  1995;  Espinosa  et  al.,  1997a).  Thus,  the  classical 
equation  of  motion  governing  the  propagation  of  a  longitudinal  wave  can  not  describe  the 
propagation  of  a  failure  wave  with  the  decrease  in  shear  strength  only,  if  other  field 
equations  are  not  invoked. 

Based  on  a  detailed  review  of  experimental,  analytical  and  computational  issues 
related  to  the  failure  wave  phenomenon,  attempts  (Chen  and  Xin,  1999;  Feng,  2000)  have 
been  made  recently  to  answer  the  fundamental  question:  what  is  the  physical  mechanism 
behind  the  failure  wave  phenomenon? 

As  indicated  in  the  previous  research,  the  jumps  of  certain  kinematic  field 
variables  involved  in  a  complete  failure  process  can  be  identified  based  on  the  transition 
between  different  governing  differential  equations  (Chen,  1996;  Chen  and  Sulsky,  1995). 
By  taking  the  initial  point  of  material  failure  as  that  point  where  the  type  of  the  governing 
differential  equations  changes,  i.e.,  a  hyperbolic  to  an  elliptic  type  for  dynamic  problems 
and  an  elliptic  to  another  elliptic  type  for  static  problems,  a  moving  material  surface  of 
discontinuity  can  be  defined  through  the  jump  forms  of  conservation  laws  across  the 
surface.  Jumps  in  density,  velocity,  strain  and  stress  can  be  accommodated  on  this 
moving  surface  of  discontinuity  between  two  material  domains.  Interestingly,  the 
problems  involving  the  type  change  of  the  governing  differential  equations,  accompanied 
by  certain  jumps  in  field  variables,  also  occur  in  other  areas  such  as  fluid  mechanics 
(Chen  and  Clark,  1995)  and  thermal  shock  wave  propagation  (Tzou,  1989  and  1997). 

Mathematically  speaking,  three  basic  kinds  of  governing  differential  equations, 
namely,  hyperbolic  (wave),  parabolic  (diffusion)  and  elliptic  (instantaneous  response 
through  the  problem  domain)  equations,  have  been  used  to  describe  different  physical 
phenomena.  However,  it  is  still  an  unsolved  challenging  problem  how  to  represent  the 
transition  among  these  three  different  kinds  of  equations  within  a  rigorous  mathematical 
framework.  If  the  transition  from  a  hyperbolic  one  to  an  elliptic  one  is  represented  via  a 
parabolic  one,  an  analytical  solution  has  been  obtained  for  a  dynamic  softening  bar  with 
the  use  of  a  simple  elastoplasticity  model  (Xin  and  Chen,  2000).  The  use  of  the  jump 
forms  of  conservation  laws,  together  with  a  simple  elastodamage  model,  also  yields  a 
diffusing  failure  front  (Chen  and  Xin,  1999).  If  the  failure  wave  were  really  a  “wave,”  its 
propagation  speed  should  increase  during  the  failure  process  because  the  material 
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stiffness  behind  the  failure  front  is  weaker  than  that  in  front  of  it  due  to  the  evolution  of 
microcracking  in  space.  As  observed  in  experiments,  however,  the  time  delay  between 
the  elastic  shock  wave  front  and  the  eruption  of  failure  appears  to  increase  with  the 
distance  into  the  material.  From  both  experimental  and  analytical  viewpoints,  therefore,  it 
appears  that  the  evolution  of  a  "failure  wave  in  space  should  be  governed  by  a  diffusion 
equation  instead  of  a  wave  equation.  Since  it  is  still  impossible  to  measure  the  real-time 
internal  failure  evolution,  it  is  difficult  to  fully  understand  the  physical  mechanism 
behind  the  failure  wave  phenomenon.  Based  on  the  experimental  data  available, 
nevertheless,  a  micromechanics-based  picture  has  been  drawn  to  explain  the  physical 
mechanism  of  failure  waves  (Feng,  2000),  as  sketched  as  follows. 

Under  plane  shock  wave  loading,  the  material  failure  below  the  HEL  occurs 
through  simultaneous  processes  of  heterogeneous  microfissuring,  shear  dilatancy  and 
void  collapsing  under  very  high  confining  stress,  which  results  in  an  increase  in  the  mean 
stress  and  a  decrease  in  the  deviatoric  stress  while  all  the  longitudinal  field  variables 
remain  unchanged.  This  particular  form  of  failure  initiates  at  the  impact  surface  where  the 
surface  defects  and  transient  loading  conditions  are  conducive  for  such  a  process,  and 
propagates  into  the  material  bulk  through  progressive  multiplication  of  microfissures  -  a 
percolation  process.  As  a  result,  the  failyre  propagation  is  of  diffusive  nature.  In  analogy 
to  the  definition  of  permanent  inelastic  strain,  the  dilated  volume  of  the  damaged  material 
at  full  release  can  be  employed  to  measure  the  extent  of  microdamage,  rendering  an 
internal  state  variable  approach  for  the  constitutive  modeling  of  failure  waves. 

A  recent  experimental  study  on  the  failure  response  of  mortar  under  impact 
loading  indicates  that  a  clearly  defined  failure  front  can  not  be  observed,  and  instead,  a 
gradual  failure  process  occurs  upon  the  arrival  of  the  loading  wave  and  propagates 
thereafter  (Grote  et  ah,  2000).  The  difference  in  the  speed  of  failure  eruption  between 
glass  and  mortar  materials  seems  to  be  related  to  the  fact  that  virgin  mortar  is 
heterogeneous  with  many  voids  as  compared  with  virgin  glass.  There  would  not  be  such  a 
difference  if  the  failure  propagation  is  not  of  diffusive  nature.  To  better  understand  the 
difference  between  different  materials,  however,  well-designed  experiments  are  still 
required  to  obtain  a  consistent  set  of  experimental  data. 

To  provide  a  general  framework  for  modeling  the  impact  response  of  certain 
brittle  solids,  a  three-dimensional  continuum  damage  model  is  developed  here,  based  on 
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the  assumption  of  shear-induced  dilatancy.  The  essential  component  of  the  proposed 
model  is  that  the  deviatoric  potential  energy  in  the  intact  material  is  converted  into  the 
volumetric  potential  energy  in  the  comminuted  and  dilated  material  during  the  time- 
dependent  failure  evolution  process.  The  progressive  percolation  of  microfissures  is 
described  by  a  nonlinear  diffusion  equation  throughout  the  continuum  body.  The 
diffusion  equation  and  wave  equation  are  then  solved  via  a  staggered  manner  in  a  single 
computational  domain.  Numerical  solutions  are  presented  and  verified  with  experimental 
data  available  to  demonstrate  the  applicability  of  the  proposed  procedure.  Future  research 
is  discussed  based  on  the  conclusion  of  the  current  work. 

2.  Constitutive  Modeling 

A  direct  notation  is  employed  here  with  boldfaced  letters  denoting  tensors  of  first 
or  higher  order.  The  Cauchy  stress  and  velocity  strain  tensors  are  used  to  implement  the 
proposed  continuum  damage  model  into  the  finite  element  program  with  an  updated 
Lagrangian  formulation,  as  discussed  later  in  the  next  section.  The  thermal  effects  are  not 
considered  here  for  the  purpose  of  simplicity. 

In  the  shocked  intact  material,  the  compressive  mean  stress  can  be  expressed  as 
om=H(n)  da) 

with 

H  =  (lb) 

V 

where  H  is  the  Hugoniot  mean  stress  response  of  the  material,  depending  on  the  volume 
compression,  /t,  in  the  elastically  shocked  state  where  the  dilated  volume  Vd  of  the 
damaged  material  is  equal  to  zero  In  Eq.  (lb),  V  and  V0  denote  the  current  and  original 
specific  volume,  respectively.  Once  the  failure  wave  is  initiated,  am  and  Vd  start  to 

increase  with  damage  at  a  fixed  value  of  fJ..  In  other  words,  there  is  no  macroscopic 
volume  change  during  the  failure  evolution,  and  the  uniaxial-strain  condition  is  preserved 
at  the  macroscopic  level  for  the  plate  impact  problem.  This  failure  response  is  modeled 
by  assuming  that 

O.  =  H(nJ  (2a) 
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with 


=  +Y±=Ysl±1l_1=P-  +  pV1  (2b) 

he  h  V  V  Po 

in  which  p  and  p0  represent  the  current  and  original  mass  density,  respectively. 
Although  Eq.  (2)  represents  a  phenomenological  approximation,  it  is  expected  to  be 
useful  as  long  as  the  damaged  material  is  sufficiently  compressed.  In  fact,  can  be 
considered  as  the  average  solid  volume  compression  of  the  damaged  material.  Based  on 
the  experimental  data  available,  the  Hugoniot  mean  stress  response  can  then  be  described 
by  a  polynomial  function  of  fie ,  namely 

® m  ~  alPe  +  a2pe  a3^e  ^ 

while  the  pressure-dependent  shear  modulus  can  be  defined  in  terms  of  a  polynomial  of 
<rm,i.e., 

G  =b0  +  bjOm  +  b2<jfn  +  bp*n  (4) 

in  which  at  and  bi  are  model  parameters  to  be  determined  from  shock  experiments.  The 
differential  form  of  Eq.  (3)  can  be  written  as 

da"‘=^LdV+^dV“  (5a) 


^  =  -(a,  +  2a2ll,+3a3f,2tjfctp-  (5b) 

=  (a,  +  2a-jit  +  3a1rf  (5c) 

With  sd=s-sv  and  ed=e-ev  denoting  deviatoric  stress  and  strain  tensors, 
respectively,  the  differential  deviatoric  stress-strain  relationship  takes  the  form  of 

dsd  =2G(ded-deid)  (6) 

in  which  ed  is  the  inelastic  deviatoric  strain  due  to  shear-induced  dilatancy.  With  i  being 
the  second  order  identity  tensor,  sv  =  -omi  and  ev  =  emi  represent  the  volumetric  parts 
of  stress  and  strain,  respectively.  The  mean  volumetric  strain,  em ,  is  defined  to  be 

3£m  =  /n~  =  3eem  +  3e\n  (7) 

vo 


6 


Imagining  that  the  material  is  compressed  from  an  initial  porous  state  (V0  +  Vd)  to  the 
current  state  (VO,  the  inelastic  mean  volumetric  strain  can  be  represented  by 

(8) 


3eL  =ln  V°-  Vd 


V0 

It  then  follows  from  Eqs.  (7)  and  (8)  that  the  elastic  mean  volumetric  strain  can  be  found 
V 


to  be  3e,„  =  In- 


V0  +  Vd 


.  Please  note  that  there  is  no  permanent  volume  compression  in  the 


kinematics  of  failure  wave  considered  here.  Damage  evolution  does  not  change  the 
current  volume,  and  instead,  changes  the  material  state  at  full  unloading.  In  other  words, 

the  condition  of  Vd  =  0  results  in  a  null  elm .  To  calculate  ed ,  assume  that  there  exists  a 


limit  surface  given  by 

f=W‘+Wd~c  =  0  (9) 

where  W*  and  Wd  are  the  volumetric  and  deviatoric  parts  of  elastic  potential  energy, 
respectively,  and  c  represents  the  limit  state  at  which  failure  occurs.  The  corresponding 
consistency  condition  can  be  written  as 

df  =  dWev  +  dWd  =  0  (10a) 


namely 

sv :  deey  +  sd  •'  de‘d  =  0  (10b) 

The  kinematics  of  failure  wave  under  plate  impact  implies  the  condition  of 
dev  =ded  =0.  In  other  words,  the  failure  wave  behind  the  shock  wave  front  will  not 
introduce  any  additional  strain  at  the  macroscopic  level,  as  discussed  before.  As  a  result, 
the  relationships  between  the  elastic  and  inelastic  parts  of  volumetric  and  deviatoric 
strains  have  the  form  of 


(11a) 


and 

deti=-deii 

respectively.  With  the  use  of  Eqs.  (10)  and  (11),  hence,  it  follows  that 


dV d  sd :  ded 


=  0 


With  the  assumption  of 


(lib) 

(12) 
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(13) 

(14) 


because  of  dVd>0  with  the  evolution  of  failure.  Thus,  the  differential  inelastic 
deviatoric  strain  tensor  takes  the  form  of 


de\  =  - — - dVd  (15) 

The  three-dimensional  isotropic  damage  model  with  shear-induced  dilatancy  is  now 


complete  for  the  given  strain  state  in  shocked  materials  under  plate  impact.  As  can  be 
seen  from  the  above  equations,  the  stress  state  depends  on  the  value  of  Vd  for  the  given 


strain  state  after  material  failure  occurs  behind  the  elastic  shock  wave. 

For  each  time  increment  At,  the  solution  steps  of  a  simple  constitutive  model 
solver  can  be  summarized  as  follows: 


2-  °m  =  alLie  +  +  a3^e  ’ 

3.  G=b0  +  blom+b2ol+bpll ; 


4.  Aeld  = 

5.  Asd  =  2G^Aed  -  Aed )  with  Aed  =  Ae  -  Aev  and  Aev  =  ev\  -  ev\  A‘ ;  and 

6.  s\  =  +  As  with  As  =  Asd+  (~Aom)i  and  Aam  =  aj  - aj  A' . 

Please  notice  that  Ae  *0  and  AVd  =  0  in  the  elastic  shock  wave  zone,  while  Ae  =  0  and 
AVd  >  0  behind  the  elastic  shock  wave  if  failure  occurs.  The  value  of  Vd  is  obtained  from 

the  damage  diffusion  equation,  as  described  later. 

To  demonstrate  the  features  of  the  proposed  damage  model,  consider  the  changes 
of  incremental  stresses  after  a  failure  wave  occurs  in  a  plate  impact  problem.  The 
incremental  longitudinal  and  lateral  stresses  are  given  by 

Asjj  =  (Asd)u  -  Aam  (16a) 

and 
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As  22  =  As  33  =  (ASj  )22  —  A(Jm  (16b) 

respectively.  Because  the  macroscopic  strain  field  is  unchanged  in  the  failure  zone  behind 
the  elastic  shock  wave,  it  follows  from  Eqs.  (5),  (6)  and  (15)  that 


Asu 


4Go_ m(sn-s22)  _£L±2£2lL±2MLAV 

3(V0+Vd)sd:sd  Vd  V 


(17a) 


=  As33=-2G(Aeid\2--^AVa 


2_9£AhiZl22lAV  ±ll£2iL±^KAV  <0 
3(V0  +  Vd)Sd:sd  d  V 


(17b) 


As  can  be  seen  from  Eq.  (17),  a  zero  change  in  the  longitudinal  stress  might  result  in  an 
increase  in  the  absolute  value  of  lateral  stresses  for  a  suitable  value  of  Vd  (Note  that 


sn  <s22  =s31<  0  in  a  plate  impact  problem).  The  effects  of  model  parameters  on  the 

failure  response  will  be  demonstrated  in  Section  4. 

The  above  damage  model  describes  only  the  evolution  of  damage  with  time  at  a 

t 

given  material  point.  The  interactions  among  different  material  points  in  a  continuum, 
which  result  in  the  evolution  of  damage  in  space,  must  be  modeled  through  an 
appropriate  manner.  Higher  order  models,  such  as  nonlocal  integral  or  strain  gradient 
models,  have  been  used  to  predict  the  evolution  of  localized  failure,  as  reviewed  by  Chen 
and  Schreyer  (1994).  However,  the  use  of  higher  order  models  yields  higher  order 
governing  differential  equations,  in  addition  to  the  ambiguity  of  additional  boundary 
conditions.  Another  approach  is  to  apply  the  jump  conditions  to  different  material 
domains  so  that  the  classical  governing  differential  equations  are  still  valid  to  simulate 
the  evolution  of  localization  (Chen  and  Sulsky,  1995;  Chen  and  Xin,  1999).  In  this  report, 
a  three-dimensional  damage  diffusion  equation  that  is  based  on  the  physical  mechanisms 
of  failure  waves  is  formulated,  as  elucidated  next,  to  predict  the  evolution  of  damage  in 
space  through  the  value  of  Vd  .  Since  both  the  wave  equation  and  diffusion  equation  are 

still  the  2nd  order  partial  differential  equations,  an  effective  numerical  procedure  can  be 
developed  with  parallel  computing  for  large-scale  model-based  simulation. 


9 


1 


As  indicated  by  Feng  (2000),  severe  loading  conditions  at  the  impact  surface, 
including  the  effect  of  impact  tilt,  may  play  the  same  role  as  the  impact  surface 
imperfection  in  initiating  the  failure  wave  in  shocked  glasses.  If  the  plane  of  entering 
wave  front  is  not  exactly  coincident  with  that  of  the  material  surface  as  it  is  usually  the 
case  in  practice,  the  transient  loading  conditions  on  the  surface  involve  a  line  of  stress 
concentration  sweeping  across  the  surface.  Because  of  the  geometry,  the  confinement 
associated  with  this  moving  stress  concentration  is  lower  than  that  of  the  uniaxial  strain 
and  the  loading  rate  can  be  significantly  higher  than  that  at  the  wave  front  inside  the 
material.  The  combination  of  the  severe  loading  conditions  and  inevitable  surface 
imperfection  may  initiate  isolated  microcracking  in  the  vicinity  of  the  impact  surface. 
Please  note  that  the  microscopic  defects,  if  any,  randomly  distributed  inside  the  specimen 
material  should  not  be  severe  enough  to  initiate  the  failure  wave,  but  will  certainly  speed 
up  the  failure  evolution  process  inside  the  material.  For  the  purpose  of  simplicity, 
however,  the  bulk  of  the  material  is  assumed  to  be  flawless  so  that  the  effect  of  randomly 
distributed  microscopic  defects  is  not  considered  here.  Based  on  the  above  analysis,  it 
appears  to  be  reasonable  to  consider  the  impact  surface  as  the  only  source  for  initiating 
the  failure  wave  in  shocked  glasses.  The  formation  of  microcracks  on  the  impact  surface 
will  in  turn  introduces  local  stress  concentrations  to  the  adjacent  downstream  material. 
Once  these  stress  concentrations  reach  a  certain  threshold,  microfissuring  will  be  initiated 
there,  and  so  on.  Since  the  limitation  of  current  computational  capabilities  does  not  allow 
a  detailed  modeling  of  the  evolution  of  each  individual  microcrack,  the  failure 
propagation  process  can  be  described  by  a  progressive  percolation  of  microfissures  in  an 
average  sense.  Based  on  the  physics  involved,  this  process  is  macroscopically  of  diffusive 
nature,  i.e.,  a  diffusion  process  starting  from  a  high  concentration  of  microcracks  to  a  low 
concentration  of  microcracks  until  a  saturated  state  is  reached  in  the  continuum  or 
discontinuous  failure  occurs.  To  activate  the  diffusion  process,  both  microscopic 
heterogeneity  and  sufficient  deviatoric  strain  energy  are  required  as  long  as  it  is  initiated. 
If  either  the  heterogeneity  or  deviatoric  strain  energy  is  not  large  enough,  the  diffusion 
process  will  not  be  active  or  die  out  eventually.  From  a  viewpoint  of  energy  conservation, 
the  evolution  of  a  failure  wave  in  space  converts  the  deviatoric  strain  energy  in  the  intact 
material  ahead  of  the  failure  front  into  the  volumetric  potential  energy  in  the  comminuted 
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material  behind  the  front.  As  a  result,  there  is  a  decrease  in  the  stress  deviator  and  an 
increase  in  the  mean  stress  in  the  uniaxial  strain  condition. 

It  should  be  pointed  out  that  even  if  the  shock  wave  compression  is  along  the 
longitudinal  direction,  the  percolation  of  microfissures  occurs  three-dimensionally. 
Therefore,  a  three-dimensional  diffusion  equation  must  be  formulated  to  govern  the 
percolation  of  microfissures,  which  is  measured  in  terms  of  the  dilated  volume,  as 
discussed  before.  In  a  standard  form,  the  diffusion  equation  in  the  three-dimensional 
space  x  with  time  t  can  be  written  as 


?lsL  =  V'[D(x,t)*VVd]  (18) 

dt 

where  D(x,t)  denotes  the  second  order  damage  diffusivity  tensor.  Although  the  rate  of 
percolation  in  the  lateral  direction  is  different  in  general  from  that  in  the  longitudinal 
direction,  the  existing  experimental  data  do  not  warrant  the  detailed  formulation  of 
D(x,t).  If  the  microscopic  details  of  percolation  in  different  orientations  are  not  pursued, 

it  is  reasonable  to  let  D(x,t)=  D(x,t)i  with 


D(x,t)  = 


0 


d 


IzIf 

yhel  ~  yf 


>0 


if  Y  <  Ythd  orVd  -0 
ifY>YTHD  or  Vd  >0 


(19) 


where  d  is  the  isotropic  diffusion  coefficient,  and 


Y  = 


:sd 


(In  a  plate  impact 


problem,  Y  =\sn  -s22 1)-  The  parameters  YHEL  and  YF  represent  the  stress  deviators  in 

the  intact  material  at  the  HEL  and  in  the  damaged  material  at  the  completion  of  failure 
process,  respectively.  As  can  be  seen  from  Eq.  (19),  the  diffusion  is  inactive  if  the  stress 
deviator  is  below  the  threshold  YTHD ,  or  if  the  dilated  volume  is  equal  to  zero.  As  long  as 
Y  >  YTHD  or  Vd  >  0 ,  the  diffusion  process  will  start  until  Y  =  YF.  Unloading  occurs 
when  Y  <YF.  Since  certain  time  is  required  for  the  lateral  microdamage  percolation  to 
finish  at  a  given  longitudinal  location,  introduce  a  time-dependent  evolution  function  into 
Eq.  (18),  which  is  defined  by 

Q[x>t)=£MY±iXdo_>o  (20) 

d  Td 
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with  Vd0  being  the  threshold  below  which  Q{x,t )  is  inactive,  and  Td  denoting  the 

characteristic  time  of  the  damage  evolution  at  the  longitudinal  location.  As  a  result,  Eq. 
(18)  then  becomes 

^  =  V  •  [D(x,t)VVd  ]+  GM  (21) 

at 

The  effect  of  lateral  percolation  on  the  longitudinal  percolation  is  reflected  through  the 
addition  of  Q{x,t). 

As  can  be  seen  from  the  above  description,  the  proposed  constitutive  model 
includes  the  equations  governing  the  strain-stress  response  at  the  local  level,  which 
mainly  consists  of  Eqs.  (2-4,  6,  15),  and  the  equations  governing  the  damage  diffusion, 
namely,  Eqs.  (19-21). 

The  propagation  of  a  failure  wave  is  simulated  through  the  evolution  in  both 
temporal  and  spatial  domains  of  the  internal  state  variable,  Vd  .  If  a  specimen  is  shocked 

beyond  the  threshold  YTHD  and  an  initial  value  of  Vd  (>Vdo)  is  assigned  to  the  vicinity 
of  the  impact  surface,  a  failure  wave  will  propagate,  behind  the  leading  elastic  shock 
wave,  through  the  specimen  that  is  assumed  to  be  originally  flawless  ( Vd  =0).  Although 

Vd  increases  during  the  damage  diffusion,  the  specific  volume  of  the  compressed 
material,  V,  remains  constant,  and  so  do  the  macroscopic  longitudinal  stress  and  particle 
velocity.  The  increase  of  Vd  results  in  the  decrease  of  the  deviatoric  potential  energy,  as 

reflected  through  the  increase  of  lateral  stresses.  In  other  words,  the  deviatoric  potential 
energy  in  the  intact  material  is  converted  into  the  volumetric  potential  energy  in  the 
comminuted  and  dilated  material  during  the  time-dependent  failure  evolution  process. 

It  is  expected  that  the  material  parameters  for  describing  the  damage  diffusion, 
namely,  d,  YTHD,  YF ,  Vd0  and  Td ,  can  be  evaluated  via  well-defined  impact  experiments 

that  measure  the  damage  evolution  profile  and  the  residual  frictional  strength  of 
comminuted  material.  To  solve  both  the  usual  equation  governing  the  shock  wave 
propagation  and  the  proposed  damage  diffusion  equation  in  a  single  computational 
domain,  a  simple  numerical  procedure  is  described  in  the  next  section. 
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3.  Numerical  Procedure 

The  numerical  procedures  for  conventional  wave  and  diffusion  equations  can  be 
found  in  the  standard  books  (Belytschko  and  Hughes,  1983;  Reddy  and  Gartling,  1994; 
among  others).  In  this  report,  central-difference  in  space  and  forward  integration  in  time 
are  used  to  solve  the  diffusion  equation,  while  constant  stress  elements  in  space  and  an 
explicit  time  integrator  are  employed  to  solve  the  wave  equation. 

To  implement  the  differential  form  of  the  proposed  damage  model  into  the  shock 
wave  code,  the  updated  Lagrangian  formulation  is  used  so  that  the  matrix  equations  of 
displacement-based  finite  elements  can  be  written  as 

Wi=W-H+«  (22) 

where  [m]  is  the  time-independent  diagonal  mass  matrix,  {«}'  the  vector  of  nodal 

accelerations  at  time  t,  {T?}'  the  vector  of  externally  applied  nodal  loads  at  t,  the 
vector  of  internal  nodal  forces  at  t  measured  with  respect  to  the  configuration  at  t,  and 


{p}‘t  the  vector  of  viscous  damping  nodal  forces  at  t  measured  with  respect  to  the 

configuration  at  t.  Thus,  the  velocity  strain  tensor  can  be  obtained  in  each  time  step  to 
find  the  corresponding  Cauchy  stress  tensor  through  the  constitutive  model.  Based  on  the 
Cauchy  stresses,  the  internal  nodal  forces  can  then  be  calculated  with  respect  to  the 
current  configuration.  To  produce  a  sharp  shock  wave  front,  the  viscous  damping  nodal 
forces  are  introduced  into  Eq.  (22),  which  are  obtained  from  element  damping  stresses. 
For  each  constant  stress  element  e  at  time  t,  the  damping  stress  takes  the  form  of 


Qe  = 


APe 


At1 


rr  Ae 

+vg  — 


Ae  |  AP  t 


h 


PePt' 


(23) 


if 


Pe~P°e 


>10  7  kg/mJ .  In  Eq.  (23),  v,  and  v  are  the  linear  and  quadratic 


coefficients  of  artificial  viscous  stress,  with  =  0.95  and  vq  =  1.0  being  chosen  for 
numerical  demonstration  in  the  next  section.  C[  is  the  current  longitudinal  wave  speed. 
The  incremental  mass  density,  the  average  mass  density  and  the  average  area  are  defined 
by 


13 


(24) 


Ape  =  p'e-p‘e  At‘ 

pr = (|  p‘, + prv  j 

4 j 

The  time  step  is  estimated  as  follows: 

=  m/njo.PAfy,  max{l.2Atl,  0.035 At  ^  (25.1) 

with 

At  i  =  At(l  -3vt^^Jvf  + 1  ~V[  (25.2) 

in  which  At  is  the  estimated  time  step  without  including  the  artificial  viscosity  in  Eq. 
(22),  and  v,  =  0.1  is  chosen  here. 

Although  the  wave  and  diffusion  equations  are  coupled  in  terms  of  Vd  after 

failure  occurs,  a  simple  procedure  is  adopted  here  to  solve  the  wave  and  diffusion 
equations  in  a  parallel  (staggered)  setting,  with  the  time  step  satisfying  both  stability 
conditions  in  a  single  computational  domain. 

* 

4.  Demonstration  and  Verification 

To  demonstrate  and  verify  the  proposed  model  and  solution  procedure,  the  shock 
experiments  chosen  by  Feng  (2000)  will  be  considered  here,  which  reveal  most  clearly 
the  quantitative  nature  of  the  failure  wave  phenomenon,  as  compared  with  other  existing 
data. 

It  should  be  pointed  out  that  there  is  no  complete  set  of  data  in  the  open  literature 
which  could  be  used  to  verify  the  model  parameters,  although  the  failure  wave 
phenomenon  has  been  observed  in  several  kinds  of  shocked  glasses.  In  the  recent 
experimental  characterization  of  the  failure  response  of  mortar  under  impact  loading 
(Grote  et  al.,  2000),  no  data  are  available  to  evaluate  Eqs.  (3)  and  (4)  because  the  elastic 
properties  are  assumed  to  be  constant  in  the  shock  response.  Hence,  the  experimental 
data,  which  have  been  collected  from  several  sources  in  a  way  as  consistent  as  possible 
(Feng,  2000),  are  used  here  to  verify  the  proposed  model. 


14 


Based  on  the  wave  profile  measurements  and  mean  stress  response  for  a  soda  lime 
glass,  the  original  mass  density,  longitudinal  and  shear  sound  speeds,  and  apparent  HEL 

stress  of  the  material  are  given  by  p0  =  2530kg  / m3 ,  cL  =  5828m/ s,  cs  =  3468m/ s 
and  oHEL  =  5.95GPa ,  respectively.  It  follows  from  a  Lagrangian  analysis  that  Eqs.  (3) 
and  (4)  become 


a m  =  45.36 pe  - 1 37.00 p2e  +288.30 p3e 

(26) 

G  =  30.43, pe  +  1.49a m  -0.60a2m  -0.35a3m 

(27) 

where  the  units  of  am  and  G  are  GPa.  The  stress  deviators  at  the  HEL  and  the  threshold 
for  initiating  the  failure  wave  are  YHEL  =  4.53GPa  and  YTHD  =  2.24GPa  (which 
corresponds  to  a  shock  stress  of  3  GPa),  respectively.  The  longitudinal  and  lateral  gauge 
data  on  a  shocked  K-8  glass  (Kanel  et  al.,  1991)  are  used  to  compare  with  the  model 
prediction.  Since  the  shock  stress  in  the  experiment  was  close  to  the  specimen  material’s 
HEL  of  about  8  GPa,  which  is  significantly  higher  than  that  of  the  model  material 
( a hel  =  5.95GPa),  a  normalization  procedure  is  used  here  to  minimize  the  effects  of 
inherent  material  property  differences.  It  is  assumed  that  the  failure  process  would  result 
in  a  50%  increase  in  the  lateral  stress,  at  which  Y  =  YE  so  that  D(x,t)-0  in  Eq.  (19),  in 
accordance  with  the  experimental  data.  Both  the  experimental  and  numerical  results  are 
then  normalized  with  respect  to  their  respective  lateral  stresses  ahead  of  the  failure  front. 
The  original  time  correlation  between  the  data  profiles  at  different  locations  is  unknown. 
For  clarity,  a  time  correlation  based  on  a  sound  speed  of  5828m/s  is  incorporated  in  the 
data.  The  model  parameters  related  to  damage  diffusion,  which  match  the  data,  are  found 

to  be  d  =  12m2 /s,  Td  =  7.7xl0~9 s  and  Vd0  =  2.0xl0~7 m3 /kg .  As  shown  in  Fig.  1, 

the  model  predicts  the  essential  feature  of  the  experimental  data.  The  longitudinal  stresses 
corresponding  to  the  lateral  stresses  are  shown  in  Fig.  2.  As  can  be  seen,  there  is  almost 
no  change  in  the  longitudinal  stress  during  the  failure  evolution.  It  should  be  pointed  out 
that  Feng  (2000)  used  an  ad  hoc  approach  in  1-D  modeling  and  simulation  so  that  the 
longitudinal  stress  must  be  fixed  to  find  the  corresponding  lateral  stress.  Without  fixing 
the  longitudinal  stress,  however,  the  3-D  damage  model  proposed  here  can  predict  the 
essential  features  of  the  observed  failure  wave  phenomenon.  The  effects  of  model 
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parameters,  which  control  the  diffusion  process,  are  illustrated  in  Figs.  3  and  4.  The 
profiles  of  Vd  and  corresponding  effective  stresses  with  time  are  shown  in  Figs.  5  and  6, 
which  are  consistent  with  the  model  formulation.  Due  to  the  effect  of  lateral  percolation, 
the  effective  stress  is  decreasing  with  the  increase  of  Vd  during  the  failure  propagation. 

To  simulate  the  transition  from  continuous  to  discontinuous  failure  modes  in  the 
future  work,  the  Material  Point  Method  (MPM),  which  was  recently  developed  for  those 
problems  such  as  penetration,  perforation,  metal  forming  and  cutting  (Sulsky  et  al., 
1994),  has  been  employed  here  to  solve  the  failure  wave  problem.  Figures  7-12 
demonstrate  the  numerical  solutions  obtained  via  the  MPM,  corresponding  to  Figs.  1-6. 
As  can  be  seen,  both  the  FEM  and  MPM  yield  the  same  solution  features  for  the  same 
problem,  although  the  spatial  discretization  procedure  used  in  the  FEM  is  different  from 
that  in  the  MPM.  It  appears  that  the  numerical  dispersion  of  the  MPM  is  more  obvious 
than  the  FEM,  a  further  discussion  on  which  is  beyond  the  scope  of  this  report. 

5.  Concluding  Remarks  and  Future  Research 

To  provide  a  general  framework  for  modeling  the  impact  response  of  certain 

t 

materials,  a  three-dimensional  continuum  damage  model  has  been  developed  in  this 
report,  based  on  the  assumption  of  shear-induced  dilatancy.  The  essential  component  of 
the  proposed  model  is  that  the  deviatoric  potential  energy  in  the  intact  material  is 
converted  into  the  volumetric  potential  energy  in  the  comminuted  and  dilated  material 
during  the  time-dependent  failure  evolution  process.  The  progressive  percolation  of 
microfissures  is  described  by  a  nonlinear  diffusion  equation  throughout  the  continuum 
body.  The  diffusion  equation  and  wave  equation  are  then  solved  via  a  staggered  manner 
in  a  single  computational  domain.  Numerical  solutions,  obtained  via  both  the  FEM  and 
MPM,  have  been  presented  and  verified  with  experimental  data  available.  It  appears  that 
the  essential  features  of  the  failure  wave  phenomenon,  as  observed  in  shock  experiments 
on  glasses,  can  be  predicted  by  the  proposed  approach. 

Although  a  considerable  progress  has  been  made  in  understanding  the  failure 
wave  phenomenon  since  it  was  reported  about  10  years  ago,  well-defined  experiments  are 
still  needed  to  provide  a  consistent  set  of  experimental  data  to  verify  the  constitutive 
modeling  and  solution  procedures.  A  combined  experimental,  analytical  and 
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computational  effort  is  a  necessity  for  us  to  fully  understand  the  physics  involved  in 
formation  and  propagation  of  failure  in  impact  problems. 
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Figure  1.  Comparison  of  the  model  prediction  with  experimental  data. 


Figure  2.  The  longitudinal  and  lateral  stresses  of  model  prediction 

corresponding  to  Fig.  1 . 
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Figure  5.  The  profiles  of  Vd  at  different  times. 
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Figure  6.  The  profiles  of  effective  stresses  corresponding  to  Fig.  5. 
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Figure  7.  Comparison  of  the  model  prediction  with  experimental 
data,  by  using  MPM. 
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Figure  8.  The  longitudinal  and  lateral  stresses  of  model  prediction 
corresponding  to  Fig.  7,  by  using  MPM. 
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Figure  9.  Effect  of  damage  diffusivity,  by  using  MPM. 


Figure  1 0.  Effect  of  the  characteristic  time  of  damage 
evolution,  by  using  MPM. 
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